Freezing of polydisperse hard spheres 
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The fluid - crystal equilibria of polydisperse mixtures of 
hard spheres have been studied by computer simulation of 
the solid phase and using an accurate equation of state for 
the fluid. A new scheme has been developed to evaluate the 
composition of crystalline phases in equilibrium with a given 
polydisperse fluid. Some common assumptions in theoretical 
approaches and their results are discussed on the light of the 
simulation results. Finally, no evidence of the existence of a 
terminal polydispersity in the fluid phase is found for poly- 
disperse hard spheres, the disagreement of this finding with 
previous molecular simulation results is explained in terms of 
the inherent limitations of some ways of modeling the chemi- 
cal potential as a function of the particle size. 



PACS numbers: 05.70.Fh, 64.70.Do, 82.70.Dd 

The phase behavior of polydisperse mixtures of hard 
spheres (PHS) has received some attention in recent 
years. Different theoretical approaches and simu- 

lation methods |4[j^ have been used to gain knowledge 
about the transition from a polydisperse fluid phase to 
crystal phase(s). The theoretical approaches use 
to involve drastic approximations regarding the compo- 
sition of the phases, and the results are often presented 
in form of stability diagrams both facts should be 

carefully taken into account when interpreting the re- 
sults. Bolhuis and Kofke Q have studied the fluid - solid 
equilibria by using molecular simulation methods, find- 
ing a "terminal polydispersity" in the fluid phase that 
they interpreted as the maximum polydispersity of a fluid 
which can originate a freezing transition, and related such 
a result with some experimental data. The origin of the 
"terminal polydispersity" in Ref jJJ^ will be addressed 
in this work. 

Let P(cr) be a given probability distribution function 
of particle diameters (PDFD). The distribution can be 
characterized by its moments, nik =< > /cqi where 
(To is a reference diameter. 

The thermodynamics of PHS fluids is very accurately 
described by the generalization of Salacuse and Stell Q to 
the polydisperse case of the equation of state (EOS) due 
to Boublik and Mansoori, Carnahan, Starling and Leland 
(BMCSL) In such an equation the pressure, p, can 

be written as: (3p ~ I3p{mi,m2,m^,ri) where ry is the 
packing fraction: 77 — -k N ra2,aQ / {QV) . N is the number 
of particles, V is the volume, and (3 = l/{kBT), with 
kB being the Boltzmann's constant and T the absolute 
temperature. 



The excess chemical potential, fi^x, in the fluid phase 
takes the form: (3^ex (c) = X]fc=o ^fecr'^, where the coeffi- 
cients Cfc depend on mi, rn2, and either rj or PpaQ. 

The goal of the present work is to evaluate the fluid 
- solid equilibrium for a given PDFD in the fluid phase. 
This point of view is the main difference with the calcu- 
lations of ref Q, however the statistical mechanics un- 
derlying both procedures is basically the same. 

In order to study polydisperse systems is convenient 
to make use of the semigrand (SG) ensemble [p|,p^, 
where the pressure, the total number of particles, N, 
and the chemical potential differences between the dif- 
ferent species and a reference one are fixed. For hard 
body interactions the basic thermodynamic differential 
relation reads, 

d [Np^io] = Vd {I3p) + Pti^dN - N^d iPfia) (1) 

where /io is the chemical potential of the reference 
species. The sum is done over the other components, Ni 
is the number of particles of species i and /i^o = fJ-i — fJ-o- 
Later, a continuous distribution of sizes will be used how- 
ever the discrete description is kept, for the shake of clar- 
ity in the equations: 

An imposed chemical potential distribution (ICPD) 
is used to perform the calculations, such a distribution 
should produce the required PDFD in the fluid phase. As 
a difference with the procedure in Ref here the ICPD 
will depend on the pressure. We have taken advantage 
of the accuracy of BMCSL EOS. Such an equation let 
us to link the fluid phase composition with the chemical 
potential distribution at a given pressure. 

Let Pa{(7) be the expected PDFD, for instance, a Gaus- 
sian distribution centered at cfq and with standard devia- 
tion (TqA. We can use as input the values of /3/i (a) given 
by: 



BMCSL 



{a,pp,X) (2) 



The actual PDFD of the fluid, P{a) will be practically 
identical to Pq due to the accuracy of the BMCSL EOS: 
The two contributions to the chemical potential can be 
grouped, by using an unique set of coefficients {ofc}: 



a-k 



k=0 



(3) 
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The coefficients Ofc will be functions of A and [3p. The 
strategy to evaluate the fluid - solid equilibrium under 
the conditions stated above lies on Gibbs-Duhem (or 
Clausius-Clapeyron) integration schemes lTo|-[l^. The 
procedure is sketched as follows: Starting for a given 
point [(3pq,\q) in which both phases are in equilibrium, 
we obtain a trajectory on the (/3p, A) plane that keep equi- 
librium conditions fulfilled. This can be done because the 
chemical potential differences can be written as functions 
of Pp and A through the coefficients ak.. Therefore, con- 
sidering N fixed: 



V 



dPp 



dm- 



dX 



(4) 



In the limit of a continuous distribution of sizes we can 
write a Clausius Clapeyron analogue equation for the co- 
existence {Pp, A) line: 



\ dX Aoe. ^ Av- ELo (dak/d i(3p))^ Amu 



(5) 



where A represents the difference between the values of 
the corresponding property in the two phases and v = 
V/N . The values of the derivatives of with respect to 
(3p and A can be evaluated numerically. 

The starting point in the Clausius Clapeyron integra- 
tion (CCI) was the monodisperse hard sphere system 
(A = 0) where the equilibrium pressure is known [ p^p4) 
to be ppa^ ~ 11.71. A second order predictor correc- 
tor has been used to advance in the integration. The 
integration step was 5X = 0.0025, the initial slope was 
found to be zero. The fluid properties have been directly 
extracted from the BMCSL EOS. A number of tests for 
several points on the (/3p, A) trajectory were performed 
by carrying out SG Monte Carlo (SGMC) simulations on 
the fluid phase using N = 256 and, within numerical ac- 
curacy, no differences between simulation and theoretical 
results were found. The solid phase was considered to be 
in a face centered cubic (FCC) ordering and its properties 
were evaluated by SGMC simulation. 

Details of the simulation procedure will be published 
elsewhere . It suffices to say that three kind of moves 
were performed, i) translation of the spheres (following 
the standard procedures), ii) changes of a particle di- 
ameter, by choosing the new diameter with probability 
proportional to exp[/3/i(a)] with u £ [0, f7,„aa;]; <^max de- 
pends on the hard sphere interactions and iii) Changes 
of volume, where we found convenient to scale simulta- 
neously the size of the particles to enhance convergence 
on the sampling. 



CCI were performed by systems with with N = 108 
and N = 256. No significant finite size effects were found 
regarding the main conclusions of the work. Some results 
for N = 256 are presented in figures. In Fig |l|the results 
of the pressure as a function of the polydispersity in the 
two phases in equilibrium. In figure |2| we plot the pack- 
ing fraction of the coexisting phases as a function of the 
polydispersity of the fluid phase. In figure |^ we show how 
the average diameter of the crystal phase increases with 
the polydispersity of the fluid phase. 

In some theoretical work |l|,^ some conjectures have 
been made regarding the possibility of finding a fluid 
in equilibrium with two or more solid phases. This re- 
sult, based on diagrams of phase stability, is not consis- 
tent with the phase rule, except in a number of singular 
points. The origin of such results lies in the theoreti- 
cal approximation that the fluid composition (polydis- 
persity) is equal to the overall composition in the solid 
phases(s). As can be seen in figures |] and |3| such condi- 
tion is not fulfilled except for very low polydispersities. 
However, one should not neglect the possibility of find- 
ing two (or more) crystalline phases which could enter 
into competition to become the solid phase in equilib- 
rium with the fluid. The phase diagram of the binary 
mixtures of hard spheres [l^ can be used as an example; 
In general for an equimolar fluid mixture x = 1/2, the 
solid phase in equilibrium with the fluid is richer in the 
large component. As the size difference increases an eu- 
tectic point appears in the phase diagrams 0, it could 
happen that for given size difference the eutectic com- 
position could become Xs = 1/2, in that case we could 
have for larger size differences, that the stable solid phase 
could become composed mainly by the small spheres. In 
order to check such a possibility for the polydisperse sys- 
tem we have made some control simulations starting from 
FCC phases composed with particles smaller that (Tq at 
pressures closed to the ones obtained in CCI. Those sys- 
tems evolved either to produce the same solid appearing 
in the CCI or to the melting of the sample. These results 
seem to discard the change of stable phase in the proce- 
dure of increasing A (at least in the range studied in this 
work). It is clear, however, that such a competition be- 
tween solid phases could appear as the freezing proceeds, 
the change of the composition of the fluid phase will alle- 
viate the phase rule restrictions, and it is quite likely to 
happen for high pressures where the fluid could even dis- 
appear as an equilibrium phase. Other possibilities have 
not been considered here, for instance, the formation of 
crystal phases with a bimodal distribution of sizes, which 
could be accommodated, for instance, in a body centered 
cubic lattice. 

Here we will discuss briefly terminal polydispersity in 
the fluid phase which appears in the results of Ref. [Q. 
In that work the ICPD has the form: 
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With this function a CCI scheme was performed from 
the monodisperse hmit (A — 0), and it was fomid that 
TOi — > when A increases and that the values of the 
reduced polydispersity, S2 — \/ m^jrn^ — 1, in the fluid 
at equihbrium with a sohd phase have to be less than 
a certain "terminal polydispersity", which the authors 
identified with some experimental results of crystalliza- 
tion of colloidal mixtures. In our simulations no such 
a terminal polydispersity appeared, however the results 
are roughly similar to those of Ref until that point. 
This apparent anomaly can be explained in terms of the 
form of ICPD. In our scheme, the excess contribution to 
the chemical potential produces a positive value of the 
coefficient 03 of the ICPD, favoring large diameters to 
compensate the effect of hard sphere repulsions, however 
the form of the ICPD given by Eq (||) produces a limit 
in the maximum packing fraction of a fluid phase with 
a certain polydispersity well below close packing. It is 
possible to estimate such a limit by computer simula- 
tion [Q or using the BMCSL EOS. In fact the terminal 
points in the diagrams of Q correspond just to the cross- 
ing between the fluid branch of the phase diagram and 
the line of such a maximum packing as a function of the 
polydispersity. As pointed out in Ref 0| the end of the 
coexistence curve is conditioned by the ICPD, however 
such an end does not seem to correspond to any relevant 
physical situation. 

The existence of a terminal polydispersity in the fluid 
phase has been also treated theoretically (See for instance 
Ref. ^), however the results are strongly influenced by 
the restrictions to the size fractionation. 

From the inspection of Fig. |l| one could think that the 
stability range of a polydisperse crystal with respect to 
the fluid can lie in a range between two pressures. This is 
not the case, we must emphasize that two points with the 
same value of the reduced polydispersity on the crystal 
"branch" do not correspond to the same PDFD (See Fig. 
^ and Fig.^) . The one corresponding to the higher pres- 
sure is associated with a distribution with greater value 
of TOi and a more negative value of the skewness p^ ] 
(i.e. the distribution has an asymmetric tail extending 
out towards small values of cr). 
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FIG. 1. Reduced pressure, Ppo'^ versus the polydispersity 
of the two phases in equilibrium. The polydispersity, S2 is 
defined as: S2 = [m2/m\ — l)^^^. For the fluid phase S2 = A. 
Continuous line represents the result for the fluid phase, sym- 
bols represent the conditions of the crystal phase at equilib- 
rium with the fluid as stated in the text 
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FIG. 2. Packing fraction of the fluid and crystal phases as 
a function of A (see the text for details) . Line and symbols as 
in Fig. 1 
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FIG. 3. Reduced average diameter of the samples, mi in 
different phases as a function of the fluid polydispersity A. 
Line and symbols as in Fig. 1 
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